SYSTEM AND METHOD OF IMAGE RECONSTRUCTION FOR 
OPTICAL TOMOGRAPHY WITH LIMITED DATA 



BACKGROUND OF THE DISCLOSURE 

1. ) Field of the Invention 

This invention relates to tomography and, more particularly, to 
a method and concomitant system wherein an image of an object is directly 
reconstructed from measurements of scattered radiation using a limited set 
of data for reconstruction. 

2. ) Description of the Background Art 

The inventive subject matter addresses the physical princi- 
ples and the associated mathematical formulations underlying the direct re- 
construction method for optical imaging in the multiple scattering regime. 
Methodologies for the direct solution to the image reconstruction problem re- 
sult. Moreover, the methodologies are generally applicable to imaging with 
any scalar wave in the diffusive multiple scattering regime and are not lim- 
ited to optical imaging. However, for the sake of elucidating the significant 
ramifications of the present inventive subject matter, it is most instructive to 
select one area of application of the methodologies so as to insure a measure 
of definiteness and concreteness to the description. Accordingly, since many 
biological systems meet the physical requirements for the application of the 
principles of the present invention, especially photon diffusion imaging prin- 
ciples, the fundamental aspects of the present inventive subject matter are 
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1 conveyed using medical imaging as an illustrative application of the method- 
ologies. 

There have been three major commericial developments in 
medical imaging that have aided in the diagnosis and treatment of numer- 
5 ous medical conditions, particularly as applied to the human anatomy; these 
developments are: (1) the Computer- Assisted Tomography (CAT) scan; (2) 
the Magnetic Resonance Imaging (MRI); and (3) the Positron Emission To- 
mography (PET) scan. 

With a CAT scanner, X-rays are transmitted through, for ex- 

10 ample, a human brain, and a computer uses X-rays detected external to the 
human head to create and display a series of images-basically cross-sections 
of the human brain. What is being imaged is the X-ray absorption function 
for unscattered, hard X-rays within the brain. CAT scans can detect, for 
instance, strokes, tumors, and cancers. With an MRI device, a computer 

15 processes data from radio signals impinging on the brain to assemble life- 
like, three-dimensional images. As with a CAT scan, such malformations as 
tumors, blood clots, and atrophied regions can be detected. With a PET 
scanner, the positions of an injected radioactive substance are detected and 
imaged as the brain uses the substance. What is being imaged is the gamma 

20 ray source position. Each of these medical imaging techniques has proved 
invaluable to the detection and diagnosing of many abnormal medical con- 
ditions. However, in many respects, none of the techniques is completely 
satisfactory for the reasons indicated in the following discussion. 
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1 In establishing optimal design parameters for a medical imag- 

ing technique, the following four specifications are most important. The 
specifications are briefly presented in overview fashion before a more de- 
tailed discussion is provided; moreover, the shortcomings of each of the con- 
5 ventional techniques are also outlined. First, it would be preferable to use a 
non-ionizing source of radiation. Second, it would be advantageous to achieve 
spatial resolution on the order of a millimeter to facilitate diagnosis. Third, 
it would be desirable to obtain metabolic information. And, fourth, it would 
be beneficial to produce imaging information in essentially real-time (on the 

10 order of one millisecond) so that moving picture-like images could be viewed. 
None of the three conventional imaging techniques is capable of achieving 
all four specifications at once. For instance, a CAT scanner is capable of 
high resolution, but it uses ionizing radiation, it is not capable of metabolic 
imaging, and its spatial resolution is borderline acceptable. Also, while MRI 

15 does use non-ionizing radiation and has acceptable resolution, MRI does not 
provide metabolic information and is not particularly fast. Finally, a PET 
scanner does provide metabolic information, but PET uses ionizing radia- 
tion, is slow, and spatial resolution is also borderline acceptable. Moreover, 
the PET technique is invasive due to the injected substance. 

20 The four specifications are now considered in more detail. With 

respect to ionizing radiation, a good deal of controversy as to its effects on the 
human body presently exists in the medical community. To ensure that the 
radiation levels are within what are now believed to be acceptable limits, PET 
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1 scans cannot be performed at close time intervals (oftentimes, it is necessary 
to wait at least 6 months between scans), and the dosage must be regulated. 
Moreover, PET is still a research tool because a cyclotron is needed to make 
the positron-emitting isotopes. Regarding spatial resolution, it is somewhat 
5 self-evident that diagnosis will be difficult without the necessary granularity 
to differentiate different structures as well as undesired conditions such as 
blood clots or tumors. With regard to metabolic information, it would be 
desirable, for example, to make a spatial map of oxygen concentration in 
the human head, or a spatial map of glucose concentration in the brain. 

10 The ability to generate such maps can teach medical personnel about disease 
as well as normal functions. Unfortunately, CAT and MRI report density 
measurements-electrons in an X-ray scanner or protons in MRI-and there 
is not a great deal of contrast to ascertain metabolic information, that is, 
it is virtually impossible to distinguish one chemical (such as glucose) from 

15 another. PET scanners have the ability to obtain metabolic information, 
which suggests the reason for the recent popularity of this technique. Finally, 
imaging is accomplished only after a substantial processing time, so real-time 
imaging is virtually impossible with the conventional techniques. 

Because of the aforementioned difficulties and limitations, there 

20 had been a major effort in the last ten years to develop techniques for gen- 
erating images of the distribution of absorption and scattering coefficients 
of living tissue that satisfy the foregoing four desiderata. Accordingly, tech- 
niques using low intensity photons would be safe. The techniques should 
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1 be fast in that optical events occur within the range of 100 nanoseconds - 
with this speed, numerous measurements could be completed and averaged 
to reduce measurement noise while still achieving the one millisecond speed 
for real-time imaging. In addition, source and detector equipment for the 
5 techniques may be arranged to produce necessary measurement data for a re- 
construction procedure utilizing appropriately-selected spatial parameters to 
thereby yield the desired one millimeter spatial resolution. Finally, metabolic 
imaging with the techniques should be realizable if imaging as localized spec- 
troscopy is envisioned in the sense that each point in the image is assigned an 

10 absorption spectrum. Such an assignment may be used, for example, to make 
a map of oxygenation by measuring the absorption spectra for hemoglobin 
at two different wavelengths, namely, a first wavelength at which hemoglobin 
is saturated, and a second wavelength at which hemoglobin is de-saturated. 
The difference of the measurements can yield a hemoglobin saturation map 

15 which can, in turn, give rise to tissue oxygenation information. 

As a result of this recent effort, a number of techniques have 
been developed and patented which satisfy the four desiderata. Representa- 
tive of the technique whereby both absorption and scattering coefficients are 
directly reconstructed is the methodology, and concomitant system, reported 

20 in U.S. Patent No. 5,747,810 ("Simultaneous Absorption and Diffusion To- 
mography System and Method Using Direct Reconstruction of Scattered Ra- 
diation" ) having Schotland as the inventor (one of the inventors of the present 

inventive subject matter). In accordance with the broad aspect of the inven- 

«> 
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1 tion in '810, the object under study is irradiated by a continuous wave source 
at a given frequency and the transmitted intensity due predominantly to dif- 
fusively scattered radiation is measured at selected locations proximate to 
the object wherein the transmitted intensity is related to both the absorp- 
5 tion and diffusion coefficients by an integral operator. The absorption and 
diffusion images of the object are directly reconstructed by executing a pre- 
scribed mathematical algorithm, determined with reference to the integral 
operator, on the transmitted intensity measurements. In addition, radiation 
at different wavelengths effects imaging as localized spectroscopy. 

10 Another technique, covered in U.S. Patent No, 5,905,261 ("Imag- 

ing System and Method Using Direct Reconstruction of Scattered Radiation" 
also having Schotland as one of the inventors) discloses and claims explicit 
inversion formulas obtained from the observation that it is possible to con- 
struct the singular value decomposition of the forward scattering operator 

15 within the diffusion approximation. In accordance with the broad aspect of 
'261 for imaging an object having variable absorption and diffusion coeffi- 
cients, the object under study is irradiated and the transmission coefficient 
due predominantly to diffusively scattered radiation is measured at appro- 
priate locations proximate to the object. The transmission intensity is re- 

20 lated to the absorption and diffusion coefficients by an integral operator. An 
image representative of the object is directly reconstructed by executing a 
prescribed mathematical algorithm, determined with reference to the inte- 
gral operator, on the transmission coefficient. The algorithm further relates 
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1 the absorption and diffusion coefficients to the transmission coefficient by a 
different integral operator. 

The art is devoid of techniques for dealing with: (1) effects of 
sampling, that is, placement of the sources and receivers in order to obtain 
5 sampled data, in the reconstruction process, as well as limiting the amount 
of data being processed; and (2) paraxial reconstruction wherein a single 
source is utilized to illuminate the object, with the scattered light being 
detected by an on-axis detector plus a small number of off- axis detectors, 
and then scanning the surface of the object with the source-detectors array 
10 while varying the frequency range of the source. 

Moreover, whereas the prior art dealt with using data in a 
posterior manner to solve the fundamental integral operator relation devised 
by the prior art, the prior art is devoid of teachings and suggestion on the use 
of both sampled and limited data which are taken into account beforehand 
15 to solve the fundamental integral operator relation. 
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1 SUMMARY OF THE INVENTION 

These and other shortcomings and limtations are obviated, in 
accordance with the present invention, by a method and concomitant system 
wherein only a limited and sampled set of data is measured and used to 
5 directly reconstruct an image of an object. 

In accordance with a broad method aspect of the present inven- 
tion, a method for generating an image of an object includes: (a) irradiating 
the object with a source of radiation, (b) measuring both a sampled and lim- 
ited data set of transmitted intensities wherein said transmitted intensities 
10 are related to at least one coefficient characterizing the image by an integral 
operator, and (c) directly reconstructing the image by executing a prescribed 
mathematical algorithm, determined with reference to said integral operator, 
on said transmitted intensities. 

In accordance with yet another broad method aspect of the 
15 present invention, the measuring includes measuring both the sampled and 
limited data set of transmitted intensities in a paraxial arrangement. 

The organization and operation of this invention will be un- 
derstood from a consideration of the detailed description of the illustrative 
embodiment, which follows, when taken in conjunction with the accompany- 
20 ing drawing. 
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l BRIEF DESCRIPTION OF THE DRAWING 

FIGS. 1A, IB, and 1C are pictorial representations for the 
cases of complete data, sampled data, and both sampled and limited data, 
respectively; 

5 FIGS. ID, IE, and IF illustrate the geometries corresponding 

to the axial, paraxial, and multiaxial cases of limited data, respectively; 

FIG. 2 depicts a sample positioned between the light source 
and detector arrangements for the slab geometry; 

FIG. 3 depicts a sample positioned between the light source 
10 and detector arrangements for the cylindrical geometry; 

FIG. 4 illustrates a high-level block diagram of an illustra- 
tive embodiment of the image reconstruction system in accordance with the 
present invention; 

FIG. 5 is a high-level flow diagram of the methodology of the 
15 present invention; and 

FIG. 6 is another high-level flow diagram of the methodology of 
the present invention based upon measurements using the paraxial geometry 
for source and detectors. 

The same element appearing in more than one FIG. has the 
20 same reference numeral. 
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l DETAILED DESCRIPTION 

To place in perspective the detailed description of the present 
inventive subject matter and thereby highlight the departure from the art as 
disclosed and claimed herein, it is both instructive and informative to first 
5 gain a basic understanding of the imaging environment in which the present 
invention operates by presenting certain foundational principles pertaining 
to the subject matter in accordance with the present invention. Accord- 
ingly, the first part of the description focuses on a high-level discussion of 
the imaging systems relevant to the inventive subject matter; this approach 
10 has the advantage of introducing notation and terminology which will aid 
in elucidating the various detailed aspects of the present invention. After 
this overview, the system aspects of the present invention, as well as the 
concomitant methodology, are presented with specificity. 
OVERVIEW OF THE PRESENT INVENTION 
15 The image reconstruction problem for optical tomography is consid- 

ered in this disclosure. The effects of sampling and limited data on this 
inverse problem is elucidated and we derive an explicit inversion which is 
computationally efficient and stable in the presence of noise. 

The propagation of near-infrared light in many biological tissues is 
20 characterized by strong multiple scattering and relatively weak absorption. 
Under these conditions, the transport of light can be regarded as occurring 
by means of diffusing waves. There has been considerable recent interest 
in the use of such waves for medical imaging. The physical problem under 

10 
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1 consideration is one of recovering the optical properties of the interior of 
an inhomogeneous medium (composed of, for example, an object or sample 
embedded in a homogeneous medium) from measurements taken on its sur- 
face. A fundamental question of substantial practical importance concerns 

5 the impact of limited data on this inverse scattering problem. This question 
arises since it is often not possible to measure all the data which is necessary 
to guarantee uniqueness or stability of a solution to the inverse problem. A 
related question is how to recover the properties of the medium with a cer- 
tain spatial resolution. Hence, it is important to understand the effects of 

10 sampling of the measured data on the quality of reconstructed images. 

To further elucidate the general notions of both sampled and limited 
data, reference is made to FIGS. 1A-C. In the direct image reconstruction of 
the prior art, the presumption has been that data measurements resulted in 
so-called "complete data" , that is, ideal data continuously measured spatially 

15 on the entire measurement surface(s). The pictorial representation 110 of 
FIG. 1A depicts this notion wherein there is a single but infinite planar 
measurement surface extending in each of the two dimensions - irregular 
shape 111 has outgoing arrows to depict the infinite extent of the planar 
representation. On the other hand, closed curve 112 depicts a limited region 

20 of the infinite planar surface over which measurements might be taken, that 
is, a finite aperature or window encompassed by closed curve 112. 

In practice, data measurements by their very nature give rise to "in- 
complete data" such that data has finite precision and is spatially limited. 



25 



11 



l FIG. IB depicts case 120 wherein spatial data is "sampled", that is, data is 
taken only at specific points on the infinite surface - one such point is shown 
by reference numeral 121. Moreover, the specific arrangement of FIG. IB is 
referred to in the sequel as a "lattice" with lattice spacing "h", that is, given 
5 a measurement point on the planar surface, said point is surrounded by four 
other points each a distance h away from the given point in the horizontal 
and vertical directions for the planar surface. 

In addition, by combining the notions of the finite window of FIG. 1A 
and the sampled data of FIG. IB, one arrives at the notion of both "limited 
10 and sampled data" case 130, which has the pictorial representation of FIG. 
1C. Aperature 112 is shown to encompass a finite number of the sampled 
data points over the planar surface. 

There are many possibe ways to obtain limited data, three of which 
are shown in FIGS. 1D-F, desribed as follows. In arrangement 140 of FIG. 
15 ID, the so-called "axial" geometry is depicted. In this geometry, a single 
light source S is placed on surface A and scattered light data is collected on 
surface B which is on-axis with source A. The source is used to illuminate 
the medium between the surfaces, wherein the medium is presumed to be 
uniform in FIG. ID. The source-detector pair is then moved to the next 
20 measurement point on the surfaces, and additional data is then collected - 
the arrow above surface A depicts the movement of the source-detector pair. 

In arrangement 150 of FIG. IE, a single source S is used to illuminate 
the medium, and the scattered light data is collected by an on-axis detector 



25 



1 (Dl) and a small number of off- axis detectors (two are shown as D2 and 
D3). The entire source-detector array is then scanned over N points on the 
surfaces, as indicated by the direction arrow. Moreover, it is often necessary 
to vary the frequency of the source over a specified range, resulting in 0(N) 

5 spatial measurements, to obtain the requisite data for image recontruction. 
We refer to this method as "scanning paraxial optical tomography" which is 
discussed in much detail later. 

Arrangement 160 of FIG. IF depicts the "multi-axial" measurement 
arrangement wherein there are N harmonically related stationary point sources 

10 on surface A, and N point detectors on opposite surface B. 

The physical problem to be solved in cases covered by FIGS. 1D-F is 
that of reconstructing the optical properties of the interior encompassed by 
the surfaces from the measurements taken on surface B. The combination of 
surface A and surface B is hereinafter called a planar "slab", so the problem 

15 to be solved may be stated as reconstructing the interior of the slab. 

In the sequel, we show that the linearized form of the corresponding 
inverse scattering problem may be solved analytically by means of an explicit 
inversion formula. This result has three important consequences. First, by 
trading spatial information for frequency information, we obtain an image 

20 reconstruction algorithm that reduces the required number of spatial mea- 
surements of the scattered field from 0(N 2 ) to O(N). Second, this algorithm 
has computational complexity that scales as N log N and is stable in the 
presence of added noise. This result should be compared with the N 2 log 
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1 N scaling of the computational complexity of the complete-data problem of 
the prior art. Third, we explicitly account for the effects of sampling of the 
data, obtaining reconstructed images whose spatial resolution scales as the 
minimum separation between the sources. 
5 In addition, when the sources and detectors are placed on a square lat- 

tice with lattice spacing h, it is possible to obtain a suitably defined solution 
to the reconstruction problem in the form of an explicit inversion formula. 
An important consequence of this result is that the fundamental limit of res- 
olution in the transverse direction is the lattice spacing h. The resolution 

10 in the depth direction is determined by numerical precision and the level of 
noise in the measurements. Resolution is further controlled by the size of the 
aperature or window (which for the square lattice is definded by W=Nh) on 
which the data is taken. 

As alluded to above, in the past decade there has been a steadily grow- 

15 ing interest in the development of optical methods for biomedical imaging. 
The near-IR spectral region is of particular importance for such applications 
because of the presence of a "window of transparency" in the absorption spec- 
trum of biological tissues between 600 and 800nm. However, the propagation 
of near-IR radiation in tissue is characterized by strong multiple scattering 

20 which renders traditional imaging methods based on ray optics inapplicable. 
Instead, the propagation of electromagnetic radiation can be described by 
the radiative transport equation (hereinafter RTE) or by the diffusion equa- 
tion (hereinafter DE). In both approaches, information about the phase of 
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the electromagnetic wave is lost and the transport of light is characterized 
by the specific intensity J(r, s) at the point r flowing in the direction s. This 
description relies on a fundamental assumption, namely, that the intensity 
rather than the amplitude of the radiation field satisfies the superposition 
principle. An important consequence of this fact is that the specific intensity 
may be expressed in terms of the Green's function as follows: 

J(r, 8) = J G(r, s; r', s')S(r', s')d 3 r f d 2 s' , (1) 

where S(r,s) is an appropriate source function and we have assumed that 
the specific intensity is stationary in time. 

In a typical experiment, light is injected into an inhomogeneous medium 
by one or more optical fibers which act as point sources. Additional fibers 
are employed for collection and subsequent detection of the transmitted light. 
Thus, if a source is located at the point r s and produces a narrow collimated 
incident beam in the direction s s , and the detector measures the specific in- 
tensity at the point flowing in the direction s^, the measurable quantity 
(up to a multiplicative constant proportional to the total power of the source 
and the efficiency of the detector) is the Green's function G(r s , s 5 ; r d , s d ). 
The inverse problem of reconstructing the optical properties of the medium 
from multiple measurements of G(v s , s s ; r d , s d ) is referred to as diffusion to- 
mography (DT). 

The approach that we will consider to the above inverse problem is 
based on the fact that G may be related to the optical properties of the 

15 



medium. This dependence, which is known to be nonlinear, significantly 
complicates the inverse problem. Indeed, the Dyson equation for G can be 
written in operator form as 

G = G 0 - G 0 VG , (2) 

where Go is the Green's function for a homogeneous medium and V is the op- 
erator which describes the deviations of the optical properties of the medium 
from their background values. From the relation G = (1 + GoV)" 1 ^, it 
can be seen that G is a nonlinear functional of V. It is possible to linearize 
the inverse problem under the assumption that V is small, as is the case in 
many physical applications. The simplest approach is to use the first Born 
approximation which is given by 

G = Go - G 0 VGo . (3) 
In this case the main equation of DT can be formulated as 

$ = G 0 VGo , (4) 

where $ = Go — G is the experimentally measurable data function. Note 
that other methods of linearization can also be used, leading to an equation 
of the form (4) with a modified expression for 3> (see section 2 below). 

Since the right-hand side of (4) contains only the unperturbed Green's 
function G 0 , the properties of G 0 are of primary importance. Although the 
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functional form of Ga can be quite complicated, as in the case of the RTE 
some useful relations may be obtained from the underlying symmetry of the 
problem. Recently we have exploited the translational invariance of the un- 
perturbed medium in the slab measurement geometry within the diffusion 
approximation. Several reconstruction algorithms have been proposed and 
numerically simulated. It was shown that taking account of translational 
invariance can result in a dramatic improvement of computational perfor- 
mance. In particular, it allows the performance of image reconstructions for 
data sets with a very large number of source detector pairs, a situation in 
which numerical reconstruction methods can not be used due to their high 
computational complexity. The effects of sampling and limited data have 
been studied and it was shown that the fundamental limit of transverse res- 
olution is given by the step size of the lattice on which sources or detectors 
are placed. Therefore, a large number of source-detector pairs is required to 
achieve the highest possible spatial resolution. 

Although many of the methods discussed in the literature may appear 
to be distinctly different, they are, in fact, special cases of a general family of 
inversion formulas which are based on certain symmetries of the unperturbed 
medium. The derivation of these results is, in some sense, independent of the 
use of diffusion approximation. Indeed, the only important property of the 
Green's function G 0 which is used is translational or rotational invariance. 
In a general curvilinear system of orthogonal coordinates x u x 2j x 3 invariance 
with respect to translation of one of the coordinates, say, x u can be mathe- 
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1 matically expressed as Gq{x u x 2} x 3 \ x' v x' 2 , x' z ) = f(x x - x\\x 2 , x 3 ; x' 2 , x' 3 ) for 
some function /. Geometries in which translational invariance exists with 
respect to two coordinates are of particular interest. For example, in the slab 
measurement geometry discussed in section 3, Go is invariant with respect to 

5 translations parallel to the measurement plane; in the cylindrical measure- 
ment geometry which is discussed in section 3.4, G 0 is invariant with respect 
to rotations about and translations parallel to the cylinder axis. 

This specification is organized as follows. In section 1 we review 
Green's functions in radiative transport and diffusion theory and their re- 

10 lation to the measurable signal. In section 2 several approaches to lineariza- 
tion of the integral equations of diffusion tomography are considered. The 
slab measurement geometry is considered in section 3. Here we discuss var- 
ious particular cases, some of which have been implemented earlier, such as 
Fourier and paraxial tomography. In section 3.4 the cylindrical measurement 

15 geometry is considered. Finally, in section 4 we give calculate the kernels for 
the integral equations considered in this specification. 
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i 1 Green's functions in radiative transport and 
diffusion theory 

We assume that all sources are harmonically modulated at a frequency u 
5 in the radio-frequency range (not to be confused with the electromagnetic 
frequency). This includes the continuous- wave (cw) regime u> = 0. In this 
case all time-dependent quantities acquire the usual factor exp(— iut), and 
the RTE in the frequency domain can be written as 

10 [s • V + (im - iuj/c)} J(r, s) - fi s J A(s, s')/(r, s')d 2 s' = e(r, s) , (5) 

where fj, t = fi a + fj, s , fx a and /i s are the absorption and scattering coefficients 
and A(s, s') is the scattering kernel with the properties A(s, s') = .4(s', s) 
and/A(s,s')d 2 s' = 1. 

In radiative transport theory, it is customary to distinguish the diffuse 
and the reduced specific intensities, denoted I d and I r , respectively. The 
reduced intensity satisfies 

s • V/ r + (ii t - iu/c)I r = s(r, s) (6) 

and the boundary conditions I r = I inc at the interface where the incident 
20 radiation with specific intensity J inc enters the scattering medium. The diffuse 
intensity I d satisfies 
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[s • V + (fr - iu/c)} I d {v, s)-f, s J A(s, s')/d(r, s')d 2 s' = e r (r, s) , (7) 
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1 where the source term due to the reduced intensity, £ r (r, s), is given by 



e r (r, s) = fjL 8 J A{s, s')I r (r, s')d 2 s' , (8) 

We will assume everywhere below that the ballistic component of the inten- 
sity given by I r at the location of detectors is negligibly small and will focus 
on the diffuse component descried by equation (7). 

An inhomogeneous medium is characterized by the spatial distribution 
/i a (r) = /z a0 +<5/i a (r) and /x s (r) = /n s o+6fi s (T). Therefore, the Green's function 
for the RTE G(r, s; r', s') satisfies the Dyson equation (2) with the operator 
V given by V = <S/i a + 6fj, 8 (l — A). The operator A with matrix elements 
(rs|A|r's') = 6(r — r')A(s,s') is assumed to be position-independent. The 
unperturbed Green's function G 0 (r, s; r', s') satisfies 



1K [i • V + (^o - uj/c)] G 0 (r 5 s; r', s') - fi sQ [ A(s,s ,/ )G 0 (r ) s // ;r , ,§ , )^ // = 
15 J 

S(t - r')<5(s - §') , (9) 

where = Mao + A^o- 

The diffusion approximation is obtained by expanding Id(r, s) to first 
order in s (for nonzero modulation frequencies, an additional condition lj <S 
20 c/t must be fulfilled, where t is defined below in (17)): 

/ d (r,s) = ^(r) + i-J.s, (10) 

where 
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1 

u(r) = - c J I d (r,s)d 2 s , J(r) = J sl d (r,s)d 2 s 
Then the density of electromagnetic energy u satisfies 



(ii) 



5 



-V • [D(r)V«(r)] + (a(r) - iuj)u(r) = S(r) , 



(12) 



and the current J is given by 



J = - 



DVu 



(13) 



where the diffusion and absorption coefficients, D and a, are given by D = 
c/3[// a + // s (l - 5)] and a = cfi a , 9 = 1$' s f A(s, s')d 2 s' is the assymetry 
parameter, and the source for the diffusion equation is given by 



The Green's function for the diffusion equation does not depend on 
the directions s and s' and we will denote its matrix elements as G(r,r'). 
The Green's function for the diffuse component of specific intensity can be 
20 obtained by applying the formula (10) to the electromagnetic energy density 
u(r) = J G(r, r')5(r')dV. Here the source term for the DE must be evaluated 
using equations (14), (15). For a point unidirectional source £(r,s) = 6(r - 
r 0 )<5(s - s 0 ) one can easily find that E(r) = /i 5 0[s o ■ (r - r 0 )] exp[-^ t s 0 • 



S(r) = E(t) - - V • DQ(r) , 



(14) 




(15) 
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(r - r 0 )]<5[r - s 0 (s 0 • r) - r 0 + s 0 (s 0 • r)] and Q(r) = gE(r)s 0 , where 9 is the 
step function and we have used the condition tu <C c/t. With the additional 
assumptions that the diffusion coefficient is homogeneous in the regions close 
to sources and detectors and g < 1, the Green's function for the diffuse 
component can be written in source-detector reciprocal form as 

G(r, s; r', s') = ± (1 + r§ • V r ) (1 - f s' ■ V r >) G(r, r') , (16) 

where 

t = 3D 0 /c . (17) 

Note that the first argument of G(r, s; r', s') corresponds to the location of the 
source and the second to the detector. Interchange of the source and detector 
positions will result in a simultaneous change of sign of s and s', so that the 
above expression conserves source-detector reciprocity. If we assume that the 
source and detector positions are on the surface 5, the above equation can 
be simplified with the use of general boundary conditions for the diffusion 
equation which are of the form 

(u + £h-Vu)\ r&s = 0, (18) 

where t is the extrapolation distance and n is an outward unit normal to 
the surface at the point r. The Green's function G(r, r') must satisfy this 
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boundary condition with respect to both of its arguments, from which it 
follows that 

G(r, ft; r', s')| r , r , 6S = ± (l - £s • n) (l + ■ n') G(v, r') . (19) 

We further assume that the source and detector optical fibers are oriented 
perpendicular to the measurement surface. Then s • n = — 1 for the source 
and s' • n' = 1 for the detector. Consequently, equation (19) takes the form 

G (^i^')U^(l + j) 2 G(r,r'). (20) 

Thus, the Green's function for the RTE has been expressed in terms of the 
Green's function for the DE and the parameters i and t. Note that in the 
limit t ~+ 0 (purely absorbing boundaries), the quantity in the right-hand 
side of (20) stays finite since G(r, r ; ) goes to zero as £ 2 for r, r' 6 S in this 
limit. 

Inhomogeneities of the medium are described in the diffusion approx- 
imation by spatial fluctuations of the absorption and diffusion coefficients: 
a(r) = a 0 +6a(r) and D(r) = D 0 +6D(r). The unperturbed Green's function 
satisfies 

(- A)V 2 + a 0 - zu;)Go(r, r') = 6(r - r') , (21) 

and the full Green's function G of the Dyson equation (2) with the interaction 
operator given by 
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V = Sa(r) - V • 6D(t)V 



(22) 



2 Linearization 

In this section we discuss linearization of integral equations for the unknown 
operator V. For simplicity, we discuss the s-independent Green's function 
for the diffusion equation, G(r, r') = (r|G|r'), and use (20) to relate it to the 
measurable signal, but the same perturbative analysis applies to the Green's 
function G(r,s;r',s'). 

In the coordinate representation, the first Born approximation (3) has 



the form 



G(r s ,r d ) = G 0 (r s ,r d ) - (r s \G 0 VG 0 \r d ) , 



(23) 



where 



{r s \G 0 VGo\r d ) = J G 0 (r s , t)V(t)G 0 (t, r d )d 3 r . 



(24) 



Consequently, the data function defined by 



0(r s , r„) = M + j) [G 0 (r s ,T d ) - G(v s , r d )} 



(25) 



satisfies the linear integral equation 



24 



0(r SJ r d ) = f 1 + J j / ^o(r 5 , r) V(r)G 0 (r, v d )d\ . (26) 

Here the factor (1 + t/i) 2 is retained for the reasons discussed in section 1 
and the constant c/Att omitted. Eq. (25) is used to calculate <j) from G, while 
in equation (26) <f> must be regarded as given and V as unknown. Note that 
(26) has the same form as (4). 

In addition to the first Born approximation, there are other ways to 
obtain an equation of the type (26) in which the measurable data function 
is linear in V. The first Rytov approximation is frequently used, in which G 
is given by 



G{r s , r d ) = G 0 (r s , r d ) exp 



(rs\G 0 VGo\r d ) 



(27) 



G 0 (r s r d ) 

Eq. (27) can be brought to the form (26) by using the following definition 
for the data function: 



0(r s ,r d ) = -[l + -J G 0 (r s ,r d )ln 



G{r s ,r d ) 



(28) 



_G 0 (r 5 ,r d ) 

Here the term under the logarithm can be identified as the transmission 
coefficient. 

Another possible approach is analogous to the mean-field approxi- 
mation. The mean field approximation is obtained from the Dyson equa- 
tion (2) by fixing the position of source and making the ansatz G(r, r d ) = 
a(r 5 , r d )G 0 (r, r d ). Substituting this expression into (2) we can formally solve 
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for a(r s , r d ) and obtain: 



G(r s ,r d ) = G 0 (r s ,r d ) 



1 + 



(r s \G 0 VG 0 \r d ) 
Go(r s , r d ) 



-l 



(29) 



Eq. (29) can be brought to the form (26) by defining the data function 
according to 



*(r„ r d ) = (l + 7) 2 [G 0 (r„ v d ) - G(r s , r d )} . (30) 

Application of different formulations of perturbation theory, as dis- 
cussed above, to calculating the Green's function for the diffusion equation 
can be qualitatively illustrated. We have considered a model situation of a 
spherical inhomogeneity of radius R characterized by the diffuse wave num- 
ber k 2 = yj a 2 j D 2 embedded in an infinite medium with the diffuse wave 
number A* = y/ati/Di, where a 1>2 and D h2 are the absorption and diffusion 
coefficients in the background medium and inside the sphere, respectively. 
The Green's function can be found in this case analytically from the scalar 
wave Mie solution for imaginary wave numbers. Placing the origin at the 
center of the sphere, we obtain G{r s ,r d ) = G 0 (r s ,r d ) - (2k 1 /irD 1 )S(r s ,r d ). 
Here the dimensionless relative shadow S(r s , r d ) is given by 

c , , ^(2Z + l)a. 

S(t„ r d ) = J2 A Pi{r, ■ r d ) , (31) 
with ai being the Mie coefficients 
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where m = k 2 /k X) Pt(x) are the Legendre polynomials, Ii(x),Ki(x) are the 
modified spherical Bessel and Hankel functions of the first kind, prime de- 
notes differentiation of functions with respect to the argument in the paren- 
thesis and the the points r 5 , r d are outside of the sphere (r s , r d > R). 

Consider plotting S as a function of the ratio k 2 /ki for different source 
detector pairs. The locations of the sources and detectors are specified in a 
Cartesian reference frame (x, y, z) with the origin in the center of the sphere 
by r s = (-L/2,0,z) and r d = (L/2,0,z) where z can take different values. 
The sphere radius was chosen to be R = 0.4L. It can be seen that, in most 
cases, the mean-field approximation is superior to the first Rytov, and the 
first Rytov is, in turn, superior to Born. It should be kept in mind that 
to obtain the contrast of the absorption or diffusion coefficient, the value of 
k 2 /ki must be squared. Thus a ten-fold increase of the absorption coefficient 
inside the sphere corresponds to k 2 /ki » 3. 

3 Planar geometry 

3.1 Integral equations in the planar geometry 

The planar geometry is illustrated in FIG. 2 by slab geometry 200. The 
medium 201 to be imaged is located between two parallel measurement planes 
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210 and 211 separated by the distance L. Intensity measurements are taken 
with multiple source-detector pairs denoted "S" (220) and "D" (221). We 
denote the coordinates of the sources and detectors as r s = (x s ,p s ) and 
rd = (Xd,p d ), respectively. Here p sd = (y s4 ,z s4 ) are two-dimensional vec- 
tors parallel to the measurement planes. Without loss of generality, we as- 
sume that "transmission" measurements are performed with x s = -L/2 and 
Xd = L/2, while "reflection" measurements with x s = x d = -L/2. A point 
inside the medium will be denoted r = (x,p), where p = (y, z ) is a two 
dimensional vector parallel to the measurement planes. Further, it is as- 
sumed that the source and detector optical fibers are oriented perpendicular 
to the measurement surfaces and that their diameters are small compared 
to all other physically relevant scales. Therefore, the measured specific in- 
tensity is given, up to a multiplicative constant, by the Green's function 
G( x s,p a ,s s = yi\x d ,p d ,s d = ±x), where plus corresponds to the transmis- 
sion geometry and minus to the reflection geometry. 

In each experiment, the parameters x s , x d , s s and s d are fixed. There- 
fore, we focus on the dependence of the data on p s , p d and u. Then with 
use of one of the linearization methods discussed in section 2, we obtain the 
integral equation 

Ps, Pd) = / r(o;, p s , p d - r) V {r)d 3 r . (33) 

A few remarks concerning the above equation are necessary. First, the func- 
tion (j) is the experimentally measurable data function. To determine <j>, it 
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is necessary to know the full Green's function G{x s , p s1 s s ; x d , p d , s d ) as well 
as the unperturbed Green's function G 0 (x s , p s , s s ; x d , p d> s d ). The latter can 
be calculated analytically, or, in some cases, measured experimentally using 
a homogeneous medium. The exact expression for </> in terms of G and G 0 
depends on the linearization method used. Second, ?7(r) is a vector repre- 
senting the deviations of optical coefficients from their background values. 
Thus, 



77(f) = 



' ^a(r) N 



(for RTE) ; ?7(r) = 



<5a(r) 



8D(r) 



(for DE) . (34) 



Correspondingly, T(w, p s , p d ; r) is a vector of functions that multiply the 
respective coefficients. The specific form of T can be found if G 0 is known; 
we will give examples of such calculations in section 4. Here we define 



TY,, n , _ I ( r *>> Ps,Pd,r),r,*Mps,P d ]r)) (for RTE), 

1 l w ) Ps> Pd) r ) - S (35) 

(r a (w, P s , p d \ r), T D (u, p s , p d - r)) (for DE) . 

A fundamental property of the kernel T is its translational invariance. 
Mathematically, this means that T(u, p s , p d ; r) depends only on p s - p and 
Pd ~ P rather than on the three parameters p s , p d and p separately, so 
that the simultaneous transformation p s d p sd + a, p -+ p + a will leave 
the kernel invariant. Therefore, T(u, p s , p d ; r) can be written as the Fourier 
integral 
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tv / d?qsd 2 qd 

r K P s » Pd\ r) = J ( 2 tt) 4 qs ' qd ' ^ GXP ' ^ P ~ Ps ^ + iqrf " ~ P ^ ' 

(36) 

where k is the vector: k = (k^k^) for the RTE and k = {k u ,k d ) for the 
DE. Note that the isotropy of space requires that k depends only on the 
absolute values of the two-dimensional vectors q S)< j. 

By introducing the new variables Ap, q and p according to p d = 
p s + Ap and q s = q + p, q d = p, we find that 



Y{u, Ap, p s ]r) = J j^K{u, Ap, q; x) exp [iq - (p - pj] , (37) 

where 



^J2«(w,q + p,p;x)exp(ip. Ap) (38) 
and the integral equation (33) takes the form 



4>(v t Ap, pj = J T(u, Ap, p s ; r)ry(r)d 3 r . (39) 

Note that in equations (37)-(39) the list of formal arguments of T and <j> 
has been changed. Thus, for example, the data function <f>{u), p s ,p s + Ap) is 
replaced by <f>(u, Ap, p 3 ). 

We now discuss the sampling of data. First, we assume that the 
sources are located on a square lattice with lattice spacing h (recall FIG. IB) 
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so that p $ = h(yn y + zn z ) where n y and n z are integers. Second, the vectors 
Ap, which specify the source-detector transverse separation, are assumed to 
belong to the set S, Ap G E. One can consider the situation when £ is a 
square lattice, commensurate with the lattice of sources, with a spacing h! as 



the whole plane. Finally, the modulation frequencies belong to a finite set 
{^ = 1,2,...,^}. 

We also consider an approach in which N d different linear combina- 
tions (with complex coefficients cy) of detector outputs are directly measured, 
allowing for the possibility of a phased-array measurement scheme. In this 
case, equations (37) and (38) must be modified according to 

/cpq 
J^2K(u, i, q; x) exp [iq ■ (p - p s )} ,i = l,...,N d (40) 



and the integral equation for phased-array measurements becomes 



a special case. Another case arises when the detectors continuously occupy 




x) °ij exp(ip • Apj) , 



(41) 



<f>{^,i,Ps) = / r(u,i,p s ;r)v(r)d 3 r , i 



N d , 



(42) 



where 



(j)(u> 



. Ps) = Y, % <H W > A Pj. Ps) : i = 1, • • • 
Ap.€S 



N d . 



(43) 
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Note that the matrix dj does not need to be square; in the ease of a con- 
tinuous set S, the summation must be replaced by. integration and by a 
vector of functions d(Ap). 

3.2 Inversion formulas 

It is convenient to define a new three-dimensional variable fj, = Ap) or 
= (<jj,i), with iS the set of such /x, and rewrite (43) as 

Mx,P,) = / "ity, P s ;rMr)d 3 r . (44) 

The integral operator T defines a map between two different Hilbert spaces 
Hi and W 2 - Eq. (44) can be written in Dirac notation as 

|0 = rfo) . (45) 

The pseudoinverse solution to (45) is given by 

\v) = T + \<i>). (46) 
Here the pseudoinverse operator is given by 

r + = (r*r)" x r* = p(rr)- 1 , (47) 

where "*" denotes Hermitian conjugation and the expressions (I^T)" 1 and 
(rr*)"" 1 must be appropriately regularized. The regularized singular-value 
decomposition (SVD) of the pseudoinverse operator is given by 
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where Q(x,e) is an appropriate regularizer, e is a small regularization pa- 
rameter, and the singular functions |/ n ) and \g n ) are eigenfunctions with 
eigenvalues u\ of the operators IT* and TT, respectively: 

rn/ n > = al\f n ) , rr\ 9n ) = al\ 9n ) . (49 ) 

In addition, the following relations hold: 

nfn)=(Tn\9n) , T\ 9n ) = a\f n ) . (50) 

To obtain the SVD for the pseudoinverse operator, we first consider 
the eigenfunctions and eigenvalues of the operator IT*. Its matrix elements 
in the basis \fip s ) are given by 

(^IIT VP' S > = / ^M M M)W) exp Hq • (p. - p' s )} . (51) 

where 

pL/2 

(#i(q)l^)= K(t*,*x)ir(iM' 9 q;x)dx. (52) 

J— L/2 

From (51), it can be seen that the effective dimensionality of the eigenproblem 
can be reduced. That is, the p s -dependent part of the eigenfunctions can be 
found analytically. Indeed, the ansatz 



33 



b*P.\fvu) = ^ exp(-m • p a ){fj\C v {u)) , (53) 

where v and u are the indexes that label the eigenfunctions with fj,, v € S, 
and u is a two-dimensional vector in the first Brillouin zone (FBZ) of the 
lattice on which the sources are placed, -n/h < u y ,u z < n/h. It can be 
verified that |/„ u ) are eigenfunctions of IT* if |C„(u)) are eigenvectors of the 
matrix M(u) defined by 

M(u) = ^M 1 (u + v) , (54) 

V 

where the v's are reciprocal lattice vectors, v = (2ir/h)(n y y + n z z). We 
denote the eigenvalues of the nonnegative definite matrix M(u) by M^(u) 
and the singular values of the problem (the eigenvalues of IT*) by a^ u . Then 
the following relations hold: 



M(u)|a(u)) = M, 2 (u)|a(u)), (55) 

rr\M = , (56) 

<r uu = h~- l M v (\i) . (57) 

Note that the singular functions |/„ u ) (53) are normalized according to 
(/t/ul/iAi') = 8 uvl 8(u - u'). 

The second set of singular functions, |^ u ), can be obtained from the 
relations (50): 
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where 



P(/x, u; Z, p) = £ # u + v; x) exp(zv • p) . (59) 

V 

To obtain an inversion formula, according to (46) and (48), we need 
to evaluate the scalar product {f U u\4>)- It can be shown by direct calculation 
that 



(/,u|0) = ^-E(^( u )W(M ) u) ) (60) 
where 4>((jl, u) is the lattice Fourier transform of <f>{n, p s ) with respect to p s : 



c£(/x, u) = E <K*t, P.) exp(iu • p s ) . (61) 

Ps 

Finally, we arrive at the following inversion formula: 



"W - ?L|^i e(w)expHu -'' ) 

x £p-(^,u;r)</i|C„(u))(C„(u)|(j'>#(M',u) . (62) 
The above result can be simplified by noting the relation 

Ee ( w) |a( tf (u)l = ■ l 63 ' 
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Using the above relation, the inversion formula (62) can be rewritten as 

V{r) = h " Lz (2^ *M-™-p) E u; r)^\M-Hn)\^W, u) . (64) 

The above formula is our main result pertaining to the planar geometry. 
In the reminder of this section, we discuss the important features of this 
inversion formula. 

First, the pseudoinverse solution (64) was derived under the assump- 
tion that the sources occupy an infinite lattice. In practice, however, the 
sources must be restricted to a finite spatial window, that is, the data is spa- 
tially limited. In this case, the inversion formula (64) becomes approximate. 
However, a high accuracy can be achieved if the edges of the window are far 
enough from the inhomogeneities of the medium. This is because the Green's 
functions in an absorbing medium (for both RTE and DE) exponentially de- 
cay in space, and the data function can be effectively replaced by zero when 
the source location is far enough from the medium inhomogeneities. 

Second, consider the computational complexity of the method. The 
variable u is continuous, but, in practice, must be discretized. The number 
of discrete vectors u should roughly correspond to the number of different 
sources used in the experiment. As discussed above, this is a finite number 
which we denote by N v We will refer to Ni as to the number of external 
degrees of freedom. Next, let the variable ^ run over N 2 discrete points. 
Here N 2 is the number of the "internal" degrees of freedom. A direct nu- 
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1 merical SVD problem will require diagonalizing a matrix with the size N1N2 
which has computational complexity 0((NiN 2 ) 3 ) operations. However, the 
inversion formula (64) requires only Ni diagonalizations of the matrix M(u) 
whose size is N 2 (hence the terms "external" and "internal" degrees of free- 

5 dom). The computational complexity of this procedure is 0(NiN$), which 
is N* times smaller than that for the purely numerical method. For large 
iVi, this is an enormous advantage. Note that one should add to the above 
estimate the number of operations necessary to Fourier-transform the data 
function and to sum over the variables /x, // and u in equation (64). The 

10 computational cost of the first task scales as O(NilogNi) with the use of 
the fast Fourier transform and, if log Ni < JV|, can be neglected. The second 
task requires 0(jVi7V|) operations, which is also negligibly small. 

Third, it can be seen from the inversion formula (64) that the trans- 
verse resolution is limited by the lattice step h. Therefore, it is practical 

15 to evaluate the function ^(r) = r)(x, p) only in such points p which coincide 
with the lattice of sources. In this case, the factor exp(zv- p) in the definition 
of P(//, u; x, p) (59) becomes equal to unity and, consequently, P becomes 
independent of p and can be written as P(/x, u;x). Then the double sum 
in (64) is a function of u and x only. This fact significantly improves the 

20 computational performance of the algorithm. 

Fourth, we discuss the case of band-limited unknown functions. So 
far we placed no restrictions on r)(x,p) except for square integrability. In 
some cases it is known a-priori that r)(x,p) is smooth, i.e., does not change 
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very fast in space. In particular, consider the case when it is known that 
the Fourier transform f}(x, q) of the function 77 turns to zero if \q y \ > ir/h or 
\q z \ > 7r/h: 

fj(x, q) = J r/(x, p) exp(zq ■ p)d 2 p = 0, if q £ FBZ . (65) 

Thus, functions which satisfy (65) are transversely band-limited to the FBZ 
of the lattice of sources. The operator T that acts on the Hilbert space of 
such functions, to the Hilbert space of the data, H 2 , can be written as 

T(/i, p s) r) = 1^ j^K{p, q; x) exp [in . (p - p s )] . (66) 

Note that integration in (66) over d 2 u is limited to the FBZ, in contrast to 
the analogous equation (37) where integration over d 2 q is carried out over 
the entire space. However, the two operators (66) and (37) are equivalent 
if they act on the space H^K This fact can be used to show that the SVD 
pseudo-inverse solution on the space of transversely band-limited functions 
7] has the form 

d 2/ LL 

V(r) = h ' /fbz (2^ GXPHU " p) £ K * {( *> U; ^^(u^W, u) , 

(67) 

where M 1 is given by (52). Thus, summation over the reciprocal lattice 
vectors that is required for calculation of P and M in the inversion formula 
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(64) is completely avoided if it is known a priori that rj is transversely band- 
limited. 

Finally, we consider the mathematical limit h — > 0, which corresponds 
to measuring the data continuously. In this case all the reciprocal lattice 
vectors become infinite, except for v = 0. Since the functions K(fi,q;x) 
decay exponentially with |q|, we have in this limit M(u) = Mi(u) and 
P(/x, u; x, p) = K(fx, u;x). We also use the relation \im h ^ 0 (h 2 (f)) = where 
<j) is the continuous Fourier transform of the data function defined by 

<K^> u ) = / d 2 p s <f>(n, p s ) exp(zu • p s ) , (68) 
to show that, in the limit h — > 0, the inversion formula (64) becomes 

v{r) = I T^^Hq-pJE^V.qi^^l^f'^I^U'.q) • (69) 

3.3 Special cases 
3.3.1 Fourier tomography 

In this section we consider the case when both sources and detectors are lo- 
cated on square lattices. We assume that the two lattices are commensurate 
and the lattice of detectors is a subset of the lattice of sources. More specifi- 
cally, let p $ = h s (yn sy +zn sz ) and p d = h d (yn dy +zn dz ) where n sy ,n sz , n dy ,n dz 
and h d /h s are integers (h d > h s ). 
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First, consider the expression (38) for K(u, Ap, q; x). For commen- 
surate lattices and h d > h s , it can be seen that Ap is on the same lattice as 
p d . Therefore, (38) can be rewritten as 

f d?w 

K(u,Ap,q;x)= — — #c(w, q + w + v d , w + v d ; x) exp(iw • Ap) . 

(70) 

Here integration is carried out over the first Brillouin zone of the lattice with 
the step h d (FBZ(/i d )) and v d = (27r/h d )(yn y + zn 2 ) is a reciprocal vector of 
the same lattice. 

Substituting (70) into (52) and (54), we obtain the following expres- 
sion for the elements of matrix M (u): 

</*|M(u)|,0 = / fbw -—-(wwlAf (uJI^/wO exp [i (w • Ap - w' • Ap')] , 

(71) 

where 



<uw|M(u)|u/w'> = £ 53 (w.w + VilAf^u + v.W + vi) (72) 

and 



(wplM^qJIwV) = / «(w,q + P,p;a:)«V.q + P',P , ;*)«fa:- (73) 
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Here v s is a reciprocal vector for the lattice with the step h s , u is in the 
FBZ of the same lattice and (ww|M(u)|a/w') can be viewed as a Fourier 
transform of (cjAp|M(u)|a/Ap') with respect to variables Ap and Ap'. 

It can be easily verified that the inverse of the matrix M(u) is given 
in terms of the inverse of M(u) by the following formula: 

MAT^u)!//) = h\ f <Pw<Pw'{ujVf\M- l {vL)\Jv,') 

«^FBZ(/i^) 

x exp [i (w • Ap - w' • Ap')] • (74) 

At the next step we substitute (74) into the inversion formula (64) 
and obtain the following result: 



" (r) = {h - h ' )2 L { ^ exp( - iu - p) 

x(a;w|M- 1 (u)|a; , w , )0(a; , ,u + w', -w') , (75) 

where 



P(uj, w, u; r) = «(^, u + v s + w + v rf , w + v d ; x) exp(zv s • p) (76) 

v s ,v d 

and 



#*>,q s ,qd)= J2 ^.P s >^)exp[z(q s -p s + q d .p d )] . (77) 

PsPi 
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1 Note that here we use the original notation <j> = <j>{uj^ p s , p d ) where p s and 
p d are the coordinates of the source and detector, respectively. 

An important feature of the obtained inversion formula is that it 
avoids numerical evaluation of the two-dimensional integral (38). Numer- 

5 ical integration according to (38) must be performed for every value of Ap, 
q and x used in the inversion formulas, which can easily become a formidable 
computational problem. However, the functions P(o;,w, u;r) and M(u) ap- 
pearing in the inversion formula (75) are expressed directly in terms of the 
functions k. The price of this simplification is that the operator M is contin- 
10 uous (unlike the discrete matrix M) and, therefore, difficult to invert. This 
problem is, however, easily avoided by replacing the integral over d 2 wd 2 w f 
by a double sum over a finite set of discrete vectors wj, Z = 1, . . . , N w \ 

v{t) = {h s h d ) 2 —r j exp(-m • p) £ £ w, u; r) 

15 JFBZ(h s ) {2iry w y w y 

x^wlM-^^l^'w^^cj'^ + w'^w') , (78) 

The expression (78) is no longer an SVD pseudo-inverse solution obtained 
on all the available data </>(uj, p 87 p d ). Instead, it can be shown that (78) 
gives the pseudo-inverse solution on the double Fourier-transformed data 
0(cj, u + w, -w) where u continuously covers FBZ(/i 5 ) while w is discrete 
and in FBZ(/i d ). 

The shortcoming of the double Fourier method is that in order to 
obtain <^(u;, u+w, — w), the data function p s , p d ) must be experimentally 
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measured for all possible points p SJ p d , even if the number of discrete vectors 
w is small. 

In general, the number of modulation frequencies used in the double 
Fourier method is arbitrary. However, two modulation frequencies (one of 
which could be zero) is always enough to simultaneously reconstruct both ab- 
sorbing and scattering inhomogeneities. If one can assume that only absorb- 
ing inhomogeneities are present in the medium, single modulation frequency 
is sufficient, which would allow one to use the continuous- wave imaging. 

Finally, it can be seen from the inversion formulas (75), (78) that the 
transverse resolution in the reconstructed images is determined by the step 
of the finer lattice (h 8 ). The discrete vectors w can be referred to as the 
"internal" degrees of freedom, similar to Ap. The number and choice of w's 
used in the reconstruction algorithm will influence the depth resolution. 

In the limit h s , h d — > 0, inversion formula (78)takes the form 

v ^ = /7^ exp (~^'^EE' ^ *( w »q + p>p; a; )( w p|^^ 1 (q)k / p , ) 

x0(u/,q + p',-p') , (79) 

where (f>(v, q s , q d ) is the continuous Fourier transform of <f)(uj, p s , p d ) with 
respect to p s and p d . 
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3.3.2 Real-space tomography 

This imaging modality is based on direct application of the inversion formula 
(64) where the set £ is assumed to contain enough points to make the inverse 
problem sufficiently determined. As in the case of double Fourier imaging, 
one or two distinct modulation frequencies can be employed. The main 
advantage of this method is that is uses real-space measurements as the input 
data and thus utilizes all experimental measurements in the most efficient 
way However, numerical evaluation of functions K(u>, Ap, q; x) according to 
the definition (38) is complicated, especially for large values of |Ap| when 
the integral is strongly oscillatory, and can lead to loss of precision. 

3.3.3 Coaxial and paraxial tomography 

The coaxial and paraxial measurement schemes are the variations of the 
real-space method with the distinction that the number of discrete vectors 
Ap which are used in the experiment is small and all the source-detector 
transverse displacements satisfy |Ap| < L. This inequality makes the nu- 
merical evaluation of the oscillatory integral (38) much easier. The additional 
degrees of freedom which are necessary for making the inverse problem well- 
defined are obtained by considering many different modulation frequencies. 
The coaxial and paraxial tomography can be used only in the "transmission" 
geometry, when sources and detectors are placed on different planes. 

In the purely coaxial measurement geometry, only one value of Ap 
is used, namely Ap = 0. Thus, the source and detector are always on axis. 
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1 If only absorbing inhomogeneities are present, it can be seen by counting 
the degrees of freedom (two for source location plus one for modulation fre- 
quency) that the inverse problem is just enough determined in this case. 
However, there is a symmetry in the integral equations which would result 

5 in appearance of certain artifacts in the reconstructed images. Namely, the 
function AT(u;,0, q;x) defined by (38) is symmetrical in x: K(u,0,q]x) = 
K(u>, 0, q; — x). Therefore, and inhomogeneity rj(x,p) and its mirror reflec- 
tion with respect to the plane x = 0, 7/ (x, p) = r)(—x, p) would produce the 
same data function (j). In this situation, the SVD pseudo-inverse solution 

10 would approach (assuming infinite numerical precision of computations and 
an infinite set of modulation frequencies) the function 77" = (rj + r] f )/2. Thus, 
if a medium has an inhomogeneity at the point (x 0 , yo, ^o), the pseudo-inverse 
solution would show an inhomogeneity at this point and at its mirror reflec- 
tion, (— rr 0) 2/o? zo)- The problem is solved by using paraxial data (recall FIG. 

15 IE) (with 0 < Ap <C L), including summation according to (40). Small 
source-detector separations of the order of one lattice step h are sufficient to 
break the symmetry and eliminate the false images. If both absorbing and 
diffusing inhomogeneities are present, at least two detectors per source are 
required to make the problem well determined. 

20 The paraxial methods are attractive due to their experimental sim- 

plicity. Indeed, instead of independently scanning the sources and detectors 
over the measurement planes, as is required in both the double Fourier and 
the real-space methods, in the paraxial measurement scheme one only needs 
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to scan a fixed source-detector "arragement" as is illustrated in FIG. IE. 



3.3.4 Plane wave detection/illumination scheme 

An especially simple reconstruction algorithm is obtained in the case when 
for each location of the source the output of all possible detectors is added up. 
Experimentally, this can be achieved with the use of a lens to either collect 
the outgoing radiation or to illuminate the medium. Both approaches are 
mathematically identical due to the source-detector reciprocity. Obviously, 
this method can be applied only in the transmission geometry, similar to 
the paraxial and coaxial methods. The plane wave illumination scheme was 
first proposed for time-resolved diffuse tomography; here we show that this 
method is a particular case of the phased-array measurement scheme which 
is discussed in section 3.1. 

We will consider a point source which is scanned over the measurement 
plane x = -L/2 and an integrated detector at x = L/2 which measures the 
quantity / d 2 p d (j)(uj, p s , p d ) . Accordingly, summation in equation (40) over 
discrete values of Ap^ must be replaced by integration over d?Ap and the 
coefficients cy replaced by unity. This results in a simple expression for 
K(u>, q; x) which is now independent of variable i: 

K(u, q; x) = k(u, q, 0; x) . (80) 

Similarly, the kernel T defined by (40) becomes independent of i: T = 
IV,p s ;r). 
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1 It can be seen that the numerical evaluation of functions K which 

is required in both real-space and paraxial modalities is altogether avoided. 
Reconstructed images are obtained by a straightforward application of for- 
mula (64) where /x = to and, similar to coaxial and paraxial imaging schemes, 

5 multiple modulation frequencies must be employed. 



3-4 Cylindrical geometry 

The cylindrical geometry is illustrated by arrangement 300 of FIG. 3. The 
data function is measured on an infinite surface R = L/2 (surface 301), 
10 where where we have used cylindrical coordinates r = ((p,z,R). The data 
function can be written as <f) = 0(cj, (p Sj z s , <fd, z<i), where (ip $) z s ) characterize 
the location of the source and (cp di Zd) of the detector. It must satisfy the 
linear equation 

15 

r2n poo [L/2 

4>{u,<Ps,Zs,Vd,z d ) = / dip dz RdRT(uj,<p s ,z s> (p d ,z d ;ip,z,R)r}((p,z,R) . 

JO J — oo JO 

(81) 

The invariance of the unperturbed medium with respect to translations along 
the 2-axis and rotations around it requires that the kernel T has a Fourier 
expansion 



T(u, <p 8 , z s , <p d , z d ] (p,z,R)= Yl J d (*2n)* K ( u > m *' 9s ' m <*' qd > R ) 
x exp [iq 3 (z - z s ) + iqd{zd - z) + im s ((p - <p a ) + im d {<p>d ~ <p)] • 



(82) 
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1 Similar to the case of planar measurements, we introduce new variables Az, 
A<£>, g, p, m and n according z d = z s + As, ^ = </> s + Ac/? and q s = q + P-> 
q d = p } m s = m + n, m d = n. We also introduce the composite variable /x, 
which in the case of the cylindrical geometry has the form fi = (w, Ac/?, As). 

5 Then the kernel T acquires the form 



r(/i, p., z s ; V,z,R) = J2 [ TTxi K ^ 9, ™; R) exp [im(<p - <p s ) + iq(z - z s )] , 

(83) 

where 

tf(/x, q,m;R) = Y, I 7^2 K ^ m + n,q + p,n,pi R) exp [i (pAz + nA<p)] . 

(84) 

Further derivations are very similar to those performed for the geom- 
etry in section 3.2. The only difference is that the variable y which in the 
case of an infinite medium varies from -oo to oo is replaced by cp which now 
varies from 0 to 2n. To build an inversion formula, we assume that z s is 
on an infinite one-dimensional lattice with step h while <p 8 takes the values 
2 7r (j - 1)/A^, j = 1,2,. . . , jV v , and consider the eigenfunctions and eigen- 
values of the operator IT*. Omitting intermediate calculations, we adduce 
the final result (analog of equation 64 for the planar geometry): 
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Here 



x^P\fx,n,u;T)(ij,\M-\n,u)\fi')4>{fj,',n,u) . (85) 



fc=-oo v 

xexp[i(vz + N v k<p)] , (86) 

OO 

M(n, u)= Y, E M i( n + ^fc, u + w ) > ( 87 ) 

fc=-oo v 
f L/2 

(n\M 1 (m,q)\ f j!)= K(fj,,q,m;R)K*(fjf,q,m;R)MR, (88) 
J o 

</>(/i, n, u) = ^ v? S) z s ) exp [t (ny? s + uz 8 )) . (89) 

All the analysis which applies to the inversion formulas in the plane 
geometry is also applicable to equation (85). For example, the inverse matrix 
M~ l (n,u) must be appropriately regularized. The inversion formula (85) 
corresponds to the real-space method in the plane geometry. However, other 
special cases can be also considered. If either A<p or Az, or both of them, 
are on a lattice (which would require that the detectors are placed on a 
lattice which is a subset of the lattice of the sources), the double Fourier 
method discussed in section 3.3.1 can be applied. Paraxial measurement 
scheme (section 3.3.3) corresponds to the case when only a few values of 
A<p = 7r + 8(p and Az are used (where 6<p < tt and Az <C L/2) in conjunction 
with multiple modulation frequencies. Similar to the planar geometry, in the 
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purely coaxial case a symmetry is present in equation (81) with respect to 
rotation of 77 (r) by the angle tt around the 2-axis. This symmetry would 
result in appearance of artifacts in the reconstructed images. The problem 
is solved by the use of off-axis data. 

The only special case discussed in section 3.3 that can not be con- 
sidered in the cylindrical geometry is the plane wave detection/illumination 
(section 3.3.4). The reason is that one can not integrate the detector output 
over all values of A<p and Az since this will necessarily include the loca- 
tion of the source. In the planar geometry sources and detectors can be 
placed on different planes which do not intersect, which is not the case in the 
cylindrical geometry. While it is possible to achieve similar mathematical 
simplifications by integrating the source and detector outputs over angles 
while both source and detector have different ^-coordinates, namely z s ^ z d 
(physically, this corresponds to using ring rather than point-like sources and 
detectors), this will lead to a complete loss of angular resolution in the recon- 
structed images. Analogously, integration of the source and detector outputs 
along infinite lines parallel to the cylinder axis and characterized by different 
angles <p 8 ^ cp d will result in a complete loss of ^-resolution. 



4 Examples of calculating the kernels 

In this section we present an example of calculating the kernel T of the in- 
tegral equations (33) and (81) and the related functions k appearing in the 
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1 Fourier expansions of T (36) and (82). We consider the diffusion approxima- 
tion in the planar and cylindrical geometries with boundary conditions given 
by (18). As discussed in section 2, regardless of the linearization method 
used, the linearized integral equation that couples the unknown operator V 

5 to the measurable data function is given by (26). Thus, to obtain expres- 
sions for k, we must calculate the unperturbed Green's functions for the DE 
in appropriate coordinates. 

4.1 Planar geometry 

In the planar geometry, the unperturbed Green's function can be written as 



Substituting this expression into the integral equation (26), where the opera- 
tor V is defined by (22), and comparing to the definition of functions k (36), 
we find that k = n D ) is expressed in terms of the functions g as 




(90) 




(91) 




1 + 1) ^ s ' qrf5 ( qs ' Xs > x )s(<hi\ x, x d ) 

dg(q s ;x s> x) dg(q d ;x,x d ) 
dx dx 



(92) 



Here an implicit dependence of the functions g on the modulation frequency 
u) is implied. 

Thus, it is sufficient to find the functions g which satisfy the DE (21) 
and the boundary conditions (18). Substituting (90) into (21), we find that 
g^x.x 1 ) must satisfy the one-dimensional equation 



, S(x - x f ) 
g(q]X,x) = — ^— L ? 



(93) 



where 



Q(q) = + fc 2 )V2 (94) 

and the diffuse wave number k is given by k 2 = (a 0 - iu)/D 0 . 

As follows from (93), the function g is a linear combination of ex- 
ponentials exp(±Qx) with coefficients depending on x' . It is continuous at 
x = x' but its first derivative experiences a discontinuity at this point: 



g(q; x! + 0, x') - g(q; x! - 0, x') = 0 , (95) 
g'(q; x' + 0, x') - g'(q; x' - 0, x') = -1/D 0 . (96) 

In addition, the boundary conditions (18) at x = ±L/2 read 



g(q;-L/2,x')-eg'(q;0,x') = 0 , 
g(q;L/2,x') + eg'(q;L/2,x') = 0 , 
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(97) 
(98) 



1 where the prime denotes differentiation with respect to x. The conditions 
(95-98) lead to the following expression for g: 



10 



1 + [Qtf] cosh [Q (L - \x - x f \)} - [l - (Ql) 2 ] 



2D 0 Q [sinh (QL) + 2Q£cosh (QL) + (Q^) 2 sinh (QL)] 
x cosh [Q (x + x f )} + 2Q£sinh [Q (L—\x — xf\)]. (99) 



This expression can be simplified if we take into account that in the integral 
equation (26) one of the arguments of the Green's functions (r or r') must 
be on the boundary. Thus, it is enough to consider the above expression in 
the limit when either x = ±L/2 or x f = ±L/2. It can be seen that these two 
limits are given by the same expression, namely 
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g{ci;x,xX =±L/2 = g(<i]x,x f )\ xt=±L/2 = —g b ( q] x,x') , (100) 



where 



ql(q' x x') ~ Sinh [Q {L ~ IX ~ gQl + Qt COSh [Q {L ~ lx - gQl finn 
' ' sinh (QL) + 2Q£ cosh (QL) + (Q£) 2 sinh (QL) ' 1 ; 

Now the functions k can be expressed in terms of the functions as 



20 «a(w,q 8 ,q d ;rr) = I 1 9b(q s ;x a ,x)g b (q d ;x,x d ) , (102) 

«i>(w,q»,qd;aO = hs-q<i9b(cis;x s ,x)g b (q d ;x,x d ) 

dx dx ' 
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The above expressions are well denned in the limits £ -> 0 and I -* oo. For 
example, for purely absorbing boundaries (t = 0) and in the transmission 
geometry (x s = -L/2 and x d = L/2), we obtain 

«a(5, q 8) q,; x)=(£X sinh [g(q s )(^/2-x)]sinh[Q(q d )(L/2 + 3 :)l 
\D 0 J sinh [Q(q s )L] sinh [Q( 

<fa)Ml04) 



MS, q s ,qd;a:) = 

(ilV [ Q(<ls)Q(q d ) cosh [Q(q s ) (1/2 - x)} cosh [Q(q d ) (L/2 + x)} 
\ D o) [ sinh [Q(q s )L] sinh [Q(q d )L] 

| q s ■ q rf sinh [Q(q s ) (L/2 - x)} sinh [(?(q rf ) (L/2 + x)] ' 
sinh [Q(q s )L]sinh[Q(q d )L] 



In the opposite limit of purely reflecting boundaries, we obtain 



(105) 



= C0S MQ(q,) (1/2 - x)] cosh [Q(q,) (L/2 + x)l 

«W sinh [Q(q s )L] sinh [Q(q s )L] ' (1 ° 6) 



n2 



sinh [Q(q,) (L/2 - x)] sinh [Q(q d ) (L/2 + x)\ 
sinh [Q(q s )L] sinh[(5(q d )L] 
+ q s • cid cosh [Q(qJ (L/2 - x)} cosh [Q(q d ) (L/2 + x)] 
Q(q,)<5(q d ) sinh [£(q s )L] sinh [Q(q d )L] 

(107) 
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i 4.2 Cylindrical geometry 

In the cylindrical geometry we use the following expansion of the unperturbed 
Green's function 

5 

f 00 r dq 

<?o(r, r') = / TKIv exp[im(v> - if/)] exp [iq(z - *')] #(m, 9 ; i?, #) . 

m=-oo *' V Z7r J 

(108) 

Similar to the planar geometry, we can express the functions k appearing in 
(82) in terms of the functions g defined above as 

10 

K a (u, m s , q s , m d , q d ; R)=(l + j^j g(m Sj q s] L/2, R)g(m d , q d ; R, L/2) , 

(109) 
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k d (u, m s , q s , m d , q d ;R) = + j 



dg(m s , g s ; L/2, R) dg{m d , q d ; R, L/2) 
dR dR 



+ [isQd + ^|r) 9(m s , q s ] L/2, R)g(m d , q d ; R, L/2) 



(110) 
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Here we took account of the fact that for both sources and detectors, R s = 
R d = L/2. 

Upon substitution of (108) into the DE (21), we find that g(m, q; R, R') 
must satisfy the one-dimensional equation 



1 d d m2 2 
RdR R dR ~R2- Qiq) 



9 (m,q;R,R') = - 6{R ^ , (111) 
UqH 
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The solution to (108) is given by a combination of modified Bessel and Han- 
kel functions of the first kind, I m (QR) and K m (QR), and is subject to the 
following conditions: 



g(m, q; 0, R') < oo , (112) 

g(m, q- R' + 0, R') - g(m, q; R' -0,R!) = 0, (113) 

g'(m, q- R' + 0, R') - g'(m, q; R' - 0, R') = -1/D 0 R' , (114) 

g(m, q; L/2, R') + tg'{m, q- L/2, R') = 0 . (115) 

The function that satisfies the above conditions is 



1 



g(m,q;R,R') = — 



K (OR )I (OR ) K m(QL/2) + QlK'JQL/2) 

Km(QR>)Im{QR<) - _____ 

xI m (QR)I m (QR') , (ii 6 ) 

where i?> and R< are the greater and lesser of R and R'. On the measurement 
surface equation (116) becomes 



g(m, q- p, L/2) = g(m, q; L/2, R) = —g b (m, q- R) , (117) 



where 



q b (m a- R) = 2 ^Q R ) 
9b {m,q,R) LIn{QL/2) + QtlML/2 y 



(118) 



Now we can express the kernel k in terms of the simpler functions g b as 
follows: 
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«a(w, m t , q s , m d , q d - R) = f j p 6 (m s , g s ; R)g b (m d , q d ; R) , (119) 



Da 



* \ 2 



%(m„ q s ; R) dg b (m d , q d ; R) 
dR dR 



( m s m d \ 
+ [Md + J 9b(m s , q s ; R)g b (m d , q d ; R) 



(120) 



Again, the above expression is well-defined in the limits I —> 0 and 
oo. Thus, for example, for purely absorbing boundaries, we have 



(2f* N 
~b~L 



Q{gs)Q{g d )i' m M<is)R]i , m Mqd)R\ 



ImM<ls)L/2)I md [Q{qd)L/2} 



™ s m d \ I ma [Q{q s )R]I md [Q(q d )R\ 



(122) 



#2 y w%.)W«JG(*)£/2]. 

In the case of purely absorbing boundaries (^ -> oo), the analogous expres- 
sions are 



K a(u,m s ,q s ,m d ,q d ;R) = 



( 2 \ 2 imAQ(qs)R)imAQ(g d )R} 

\D 0 Lj Q(q s )Q(q d )l' ms [Q(^)L/2}I^[Q(q d )L/2) ' 

(123) 



k d (u, m s , q s ,m d , q d ; R) = 

\D 0 i 



UQ(q s WmAQ(q d )R) 
[UQ(qs)LmUQ(^)L/2] 



R'q s q d + m a m d 7 m , [Q(q s )R]I md [Q(q d )R] 

R 2 Q(q s )Q(q d ) UQ(qs)L/2]i' md [Q(q d )L/2}_ 



• (124) 
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l 5. SYSTEM DIAGRAM 

As depicted in high-level block diagram form in FIG. 4, sys- 
tem 400 is an optical tomography system for generating an image of an 
scatter er /object using scattering data measured when emanating from an 
5 object in response to a source illuminating the object. The measurements 
which are obtained result in a sampled and limited set of scattering data. 
In particular, object 201 is shown as being under investigation. System 400 
is composed of: source arrangement 420 for probing the object 201; data 
acquisition detector arrangement 430 for detecting the scattering data corre- 

10 sponding to scattering from object 201 at one or more locations proximate to 
object 201; position controller 440 for controlling the locations of detector (s) 
430 and source(s) 420; and computer processor 450, having associated input 
device 460 (e.g., a keyboard) and output device 470 (e.g., a graphical display 
terminal). Computer processor 450 has as its inputs positional information 

15 from controller 440 and the measured scattering data from detector (s) 430. 

Computer 450 stores a computer program which implements 
the direct reconstruction algorithm; in particular, the stored program pro- 
cesses the measured scattering data to produce the image of the object under 
study using a prescribed mathematical algorithm. The algorithm is, gener- 

20 ally, based upon both a sampled and limited data set of measurements, with 
the data in a preferred embodiment being obtained using the paraxial mode 
of measurement. The algorithm selected for reconstruction depends upon 
the measurement geometry and the particular source-detector arrangement. 
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Each corresponding algorithm has been described in detail above. 

6. FLOW DIAGRAM 

An embodiment illustrative of the methodology carried out by 
the subject matter of the present invention is set forth in high-level flow dia- 
gram 500 of FIG. 5 in terms of the illustrative system embodiment shown in 
FIG. 4. With reference to FIG. 5, the processing effected by block 510 enables 
source 420 and data acquisition detector 430 so as to measure the scattering 
data, which is both sampled and limited, emanating from scatterer 201 due to 
illumination from source 420. These measurements are passed to computer 
processor 450 from data acquisition detector 430 via bus 431. Next, pro- 
cessing block 520 is invoked to compute the linear operator for the sampled 
and limited data set. In turn, processing block 530 is operated to execute 
the reconstruction algorithm using the linear operator formulation. Finally, 
as depicted by processing block 540, the reconstructed image is provided to 
output device 470 in a form determined by the user; device 470 may be, for 
example, a display monitor or a more sophisticated three-dimensional display 
device. 

Another embodiment illustrative of the methodology carried 
out by the subject matter of the present invention is set forth in high-level 
flow diagram 600 of FIG. 6 in terms of the illustrative system embodiment 
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shown in FIG. 4. With reference to FIG. 6, the processing effected by block 
610 enables source 420 and data acquisition detector 430 so as to measure the 
scattering data emanating from scatterer 201 due to illumination from source 
420. The scattering data, which is both sampled and limited, is produced 
using a paraxial measurement configuration. These measurements are passed 
to computer processor 450 from data acquisition detector 430 via bus 431. 
Next, processing block 620 is invoked to compute the linear operator for the 
sampled and limited data set. In turn, processing block 630 is operated to 
execute the reconstruction algorithm using the linear operator formulation. 
Finally, as depicted by processing block 640, the reconstructed image is pro- 
vided to output device 470 in a form determined by the user; device 470 may 
be, for example, a display monitor or a more sophisticated three-dimensional 
display device. 

Although the present invention has been shown and described 
in detail herein, those skilled in the art can readily devise many other var- 
ied embodiments that still incorporate these teachings. Thus, the previous 
description merely illustrates the principles of the invention. It will thus be 
appreciated that those with ordinary skill in the art will be able to devise var- 
ious arrangements which, although not explicitly described or shown herein, 
embody principles of the invention and are included within its spirit and 
scope. Furthermore, all examples and conditional language recited herein 
are principally intended expressly to be only for pedagogical purposes to aid 
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the reader in understanding the principles of the invention and the concepts 
contributed by the inventor to furthering the art, and are to be construed as 
being without limitation to such specifically recited examples and conditions. 
Moreover, all statements herein reciting principles, aspects, and embodiments 
of the invention, as well as specific examples thereof, are intended to encom- 
pass both structural and functional equivalents thereof. Additionally, it is 
intended that such equivalents include both currently know equivalents as 
well as equivalents developed in the future, that is, any elements developed 
that perform the function, regardless of structure. 

In addition, it will be appreciated by those with ordinary skill 
in the art that the block diagrams herein represent conceptual views of illus- 
trative circuitry embodying the principles of the invention. 
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